Differential expression no well 4
Maddie Lombardo
1/14/2022
reading the functions in that i’ll use
source(here::here('run_standard_deseq.R'))
source(here::here('make_volcano_plot.R'))reading in the files
feature_counts_file_path <- file.path(here::here("data","feature_counts_bdnf_No4"))running a deseq 1hr
one_hour_dds <- run_standard_deseq(feature_counts_file_path, base_grep = "CONTROL",
contrast_grep = "BDNF",
grep_pattern = "1hr",
baseName = "control",
contrastName = 'BDNF')## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians## The following object is masked from 'package:dplyr':
##
## count##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## The following object is masked from 'package:Biobase':
##
## rowMedians##
## Attaching package: 'data.table'## The following object is masked from 'package:SummarizedExperiment':
##
## shift## The following object is masked from 'package:GenomicRanges':
##
## shift## The following object is masked from 'package:IRanges':
##
## shift## The following objects are masked from 'package:S4Vectors':
##
## first, second## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transpose## [1] "BDNF_1B_1hr" "BDNF_2B_1hr" "BDNF_3B_1hr" "CONTROL_1A_1hr"
## [5] "CONTROL_2A_1hr" "CONTROL_3A_1hr"
## [1] "BDNF_1B_1hr" "BDNF_2B_1hr" "BDNF_3B_1hr" "CONTROL_1A_1hr"
## [5] "CONTROL_2A_1hr" "CONTROL_3A_1hr"
## [1] "Filtered Genes by CPM greater than 0.5 in a least 2 samples"
## keep
## FALSE TRUE
## 35930 24739## converting counts to integer mode## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing##
## out of 24739 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 297, 1.2%
## LFC < 0 (down) : 145, 0.59%
## outliers [1] : 0, 0%
## low counts [2] : 11511, 47%
## (mean count < 58)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results## Joining, by = "Geneid"make a volcano plot 1hr
one_hour_volcano <- make_volcano_plot(one_hour_dds$results_table)## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.deseq 2hr
two_hour_dds <- run_standard_deseq(feature_counts_file_path, base_grep = "CONTROL",
contrast_grep = "BDNF",
grep_pattern = "2hr",
baseName = "control",
contrastName = 'BDNF')## [1] "BDNF_1D_2hr" "BDNF_2D_2hr" "BDNF_3D_2hr" "CONTROL_1C_2hr"
## [5] "CONTROL_2C_2hr" "CONTROL_3C_2hr"
## [1] "BDNF_1D_2hr" "BDNF_2D_2hr" "BDNF_3D_2hr" "CONTROL_1C_2hr"
## [5] "CONTROL_2C_2hr" "CONTROL_3C_2hr"
## [1] "Filtered Genes by CPM greater than 0.5 in a least 2 samples"
## keep
## FALSE TRUE
## 27648 33021## converting counts to integer mode## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing##
## out of 33021 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 1420, 4.3%
## LFC < 0 (down) : 6390, 19%
## outliers [1] : 0, 0%
## low counts [2] : 14085, 43%
## (mean count < 20)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results## Joining, by = "Geneid"volcano 2hr
two_hr_volcano <- make_volcano_plot(two_hour_dds$results_table)## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.deseq 6h
six_hr_dds <- run_standard_deseq(feature_counts_file_path, base_grep = "CONTROL",
contrast_grep = "BDNF",
grep_pattern = "6hr",
baseName = "control",
contrastName = 'BDNF')## [1] "BDNF_1F_6hr" "BDNF_2F_6hr" "BDNF_3F_6hr" "CONTROL_1E_6hr"
## [5] "CONTROL_2E_6hr" "CONTROL_3E_6hr"
## [1] "BDNF_1F_6hr" "BDNF_2F_6hr" "BDNF_3F_6hr" "CONTROL_1E_6hr"
## [5] "CONTROL_2E_6hr" "CONTROL_3E_6hr"
## [1] "Filtered Genes by CPM greater than 0.5 in a least 2 samples"
## keep
## FALSE TRUE
## 27726 32943## converting counts to integer mode## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing##
## out of 32943 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 963, 2.9%
## LFC < 0 (down) : 2487, 7.5%
## outliers [1] : 2, 0.0061%
## low counts [2] : 17245, 52%
## (mean count < 31)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results## Joining, by = "Geneid"volcano 6h
six_hr_volcano <- make_volcano_plot(six_hr_dds$results_table)## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.now we are going to use the deseq data to do gene ontology hour 1
results_1hr <- one_hour_dds$results_table
up_1hr_ensembl <- results_1hr %>%
dplyr::filter(padj < 0.05 & log2FoldChange > 0 ) %>%
mutate(ensembl = gsub("\\..*", "", Geneid)) %>%
arrange(-log2FoldChange) %>%
pull(ensembl)
down_1hr_ensembl <- results_1hr %>%
dplyr::filter(padj < 0.05 & log2FoldChange < 0 ) %>%
mutate(ensembl = gsub("\\..*", "", Geneid)) %>%
arrange(log2FoldChange) %>%
pull(ensembl)
gp_up_1hr = gost(up_1hr_ensembl, organism = "hsapiens", ordered_query = TRUE)
gostplot(gp_up_1hr)gp_down_1hr = gost(down_1hr_ensembl, organism = "hsapiens", ordered_query = TRUE)
gostplot(gp_down_1hr)GO hr 2
results_2hr <- two_hour_dds$results_table
up_2hr_ensembl <- results_2hr %>%
filter(padj < 0.05 & log2FoldChange > 0 ) %>%
mutate(ensembl = gsub("\\..*", "", Geneid)) %>%
pull(ensembl)
down_2hr_ensembl <- results_2hr %>%
filter(padj < 0.05 & log2FoldChange < 0 ) %>%
mutate(ensembl = gsub("\\..*", "", Geneid)) %>%
pull(ensembl)
gp_up_2hr <- gost(up_2hr_ensembl, organism = "hsapiens")
gp_down_2hr <- gost(down_2hr_ensembl, organism = "hsapiens")go 6h
results_6hr <- six_hr_dds$results_table
up_6hr_ensembl <- results_6hr %>%
filter(padj < 0.05 & log2FoldChange > 0 ) %>%
mutate(ensembl = gsub("\\..*", "", Geneid)) %>%
pull(ensembl)
down_6hr_ensembl <- results_6hr %>%
filter(padj < 0.05 & log2FoldChange < 0 ) %>%
mutate(ensembl = gsub("\\..*", "", Geneid)) %>%
pull(ensembl)
gp_up_6hr = gost(up_6hr_ensembl, organism = "hsapiens")
gp_down_6hr = gost(down_6hr_ensembl, organism = "hsapiens")bg_ids_1hr <- one_hour_dds$results_table$Geneid %>% gsub("\\..*", "", .)
e_go_cc_1hr_up <- enrichGO(gene = up_1hr_ensembl,
universe = bg_ids_1hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
e_go_mf_1hr_up <- enrichGO(gene = up_1hr_ensembl,
universe = bg_ids_1hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
e_go_bp_1hr_up <- enrichGO(gene = up_1hr_ensembl,
universe = bg_ids_1hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
if(any(e_go_cc_1hr_up@result$p.adjust <0.01)){
cnetplot(e_go_cc_1hr_up)
dotplot(e_go_cc_1hr_up)
}
if(any(e_go_bp_1hr_up@result$p.adjust <0.01)){
cnetplot(e_go_bp_1hr_up)
dotplot(e_go_bp_1hr_up)}
bg_ids_2hr <- two_hour_dds$results_table$Geneid %>% gsub("\\..*", "", .)
e_go_cc_2hr_up <- enrichGO(gene = up_2hr_ensembl,
universe = bg_ids_2hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
e_go_mf_2hr_up <- enrichGO(gene = up_2hr_ensembl,
universe = bg_ids_2hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
e_go_bp_2hr_up <- enrichGO(gene = up_2hr_ensembl,
universe = bg_ids_2hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
if(any(e_go_cc_2hr_up@result$p.adjust <0.01)){
cnetplot(e_go_cc_2hr_up)
dotplot(e_go_cc_2hr_up)
}
if(any(e_go_bp_2hr_up@result$p.adjust <0.01)){
cnetplot(e_go_bp_2hr_up)
dotplot(e_go_bp_2hr_up)}
bg_ids_6hr <- six_hr_dds$results_table$Geneid %>% gsub("\\..*", "", .)
e_go_cc_6hr_up <- enrichGO(gene = up_6hr_ensembl,
universe = bg_ids_6hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
e_go_mf_6hr_up <- enrichGO(gene = up_6hr_ensembl,
universe = bg_ids_6hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
e_go_bp_6hr_up <- enrichGO(gene = up_6hr_ensembl,
universe = bg_ids_6hr,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
if(any(e_go_cc_6hr_up@result$p.adjust <0.01)){
cnetplot(e_go_cc_6hr_up)
dotplot(e_go_cc_6hr_up)
}
if(any(e_go_bp_6hr_up@result$p.adjust <0.01)){
cnetplot(e_go_bp_6hr_up)
dotplot(e_go_bp_6hr_up)}